#### R script for PEC main text code
#### Jinhyuk, 04/11/2021
###########################################

###############
### Setup R ###
###############

### Clear terminal
cat("\014")

### Clear space
rm(list = ls())

### library packages

#install.packages('rio')
library(rio)
#install.packages('dplyr')
library(dplyr)
#install.packages('mlr')
library(mlr)
#install.packages('lattice')
library(lattice)
#install.packages('ggplot2')
library(ggplot2)
#install.packages('ggthemes')
library(ggthemes)
#install.packages('gridExtra')
library(gridExtra)
#install.packages('sandwich')
library(sandwich)
#install.packages('lmtest')
library(lmtest)
#install.packages('lmPerm')
library(lmPerm)
#install.packages('coefplot')
library(coefplot)
#install.packages('stargazer')
library(stargazer)
#install.packages('hrbrthemes')
library(hrbrthemes)
#install.packages('officer')
library(officer)
#install.packages('flextable')
library(flextable)
#install.packages('meta')
library(meta)
#install.packages('grep')
library(grep)
#install.packages('MASS')
library(MASS)

### Load custom functions
#### Multiplot function
# ggplot objects can be passed in ..., or to plotlist (as a list of ggplot objects)
# - cols:   Number of columns in layout
# - layout: A matrix specifying the layout. If present, 'cols' is ignored.
#
# If the layout is something like matrix(c(1,2,3,3), nrow=2, byrow=TRUE),
# then plot 1 will go in the upper left, 2 will go in the upper right, and
# 3 will go all the way across the bottom.
#
multiplot <- function(..., plotlist=NULL, file, cols=1, layout=NULL) {
  library(grid)
  
  # Make a list from the ... arguments and plotlist
  plots <- c(list(...), plotlist)
  
  numPlots = length(plots)
  
  # If layout is NULL, then use 'cols' to determine layout
  if (is.null(layout)) {
    # Make the panel
    # ncol: Number of columns of plots
    # nrow: Number of rows needed, calculated from # of cols
    layout <- matrix(seq(1, cols * ceiling(numPlots/cols)),
                     ncol = cols, nrow = ceiling(numPlots/cols))
  }
  
  if (numPlots==1) {
    print(plots[[1]])
    
  } else {
    # Set up the page
    grid.newpage()
    pushViewport(viewport(layout = grid.layout(nrow(layout), ncol(layout))))
    
    # Make each plot, in the correct location
    for (i in 1:numPlots) {
      # Get the i,j matrix positions of the regions that contain this subplot
      matchidx <- as.data.frame(which(layout == i, arr.ind = TRUE))
      
      print(plots[[i]], vp = viewport(layout.pos.row = matchidx$row,
                                      layout.pos.col = matchidx$col))
    }
  }
}

### Load data
df <- read.csv("df.csv")
df2 <- read.csv("df2.csv")

#### Impact of coalition forming on voters' perception about a party's ideological placement

### Experiment 1: Komeito

## Figure 1: Komeito with socio-demographic covariates (DV=Ideological Placement)

mod.1 <- lm(komeito_placement ~ one_noc + one_llc + one_rlc_1 + one_rlc_2 + one_rlc_3 + female + as.numeric(Q2.1) + edu + as.numeric(Q13_1) + income, data = df)
summary(mod.1)

stargazer(mod.1)

# Plot 

coef.1<-summary(mod.1)$coefficients[,1] # coef
se.1<-summary(mod.1)$coefficients[,2] # se

mod.df.1<-data.frame(cbind(coef.1, se.1))
mod.df.1<-mod.df.1[-c(1,7,8,9,10,11),]
varnames.1<-c("No Coalition","Left","Right1","Right2","Right3")

mod.df.1$upper90 <- mod.df.1$coef.1 + qnorm(0.95)*mod.df.1$se.1
mod.df.1$lower90 <- mod.df.1$coef.1 - qnorm(0.95)*mod.df.1$se.1
mod.df.1$upper95 <- mod.df.1$coef.1 + qnorm(0.975)*mod.df.1$se.1
mod.df.1$lower95 <- mod.df.1$coef.1 - qnorm(0.975)*mod.df.1$se.1
mod.df.1$varnames.1 <- varnames.1

p1.1<-ggplot() + 
  geom_linerange(data=mod.df.1, mapping=aes(x=varnames.1, ymin=lower90, ymax=upper90), size=1.5) + 
  geom_linerange(data=mod.df.1, mapping=aes(x=varnames.1, ymin=lower95, ymax=upper95), size=0.8) + 
  geom_point(data=mod.df.1, mapping=aes(x=varnames.1, y=coef.1), size=3) + 
  geom_hline(yintercept=0, colour="red", size=0.5) + 
  labs(x="Variable Names", y="Coefficient", title="Ideological Placement: Komeito (w/ Covariates)") +  # Labels
  coord_flip() +  # Rotate the plot
  theme_bw()  # Nicer theme
p1.1

### Experiment 2: DPFP

## Figure 2: DPFP with socio-demographic covariates (DV=Ideological Placement)

mod2.1 <- lm(dpfp_placement ~ two_noc + two_llc + two_rlc + female + as.numeric(Q2.2) + edu + as.numeric(Q13_1) + income, data = df2)
summary(mod2.1)

stargazer(mod2.1)

# Plot

coef2.1<-summary(mod2.1)$coefficients[,1] # coef
se2.1<-summary(mod2.1)$coefficients[,2] # se

mod.df2.1<-data.frame(cbind(coef2.1, se2.1))
mod.df2.1<-mod.df2.1[-c(1,5,6,7,8,9),]
varnames2.1<-c("No Coalition","Left","Right")

mod.df2.1$upper90 <- mod.df2.1$coef2.1 + qnorm(0.95)*mod.df2.1$se2.1
mod.df2.1$lower90 <- mod.df2.1$coef2.1 - qnorm(0.95)*mod.df2.1$se2.1
mod.df2.1$upper95 <- mod.df2.1$coef2.1 + qnorm(0.975)*mod.df2.1$se2.1
mod.df2.1$lower95 <- mod.df2.1$coef2.1 - qnorm(0.975)*mod.df2.1$se2.1
mod.df2.1$varnames2.1 <- varnames2.1

p2.1<-ggplot() + 
  geom_linerange(data=mod.df2.1, mapping=aes(x=varnames2.1, ymin=lower90, ymax=upper90), size=1.5) + 
  geom_linerange(data=mod.df2.1, mapping=aes(x=varnames2.1, ymin=lower95, ymax=upper95), size=0.8) + 
  geom_point(data=mod.df2.1, mapping=aes(x=varnames2.1, y=coef2.1), size=2) + 
  geom_hline(yintercept=0, colour="red", size=0.5) + 
  labs(x="Variable Names", y="Coefficient", title="Ideological Placement: DPFP (w/ Covariates)") +  # Labels
  coord_flip() +  # Rotate the plot
  theme_bw()  # Nicer theme
p2.1

###########################################################################################

#### Impact of coalition forming on voters' perception about the likelihood of contruction of a particular coalition

### Experiment 1: Komeito

## Figure 3

# Likelihood of no coalition formation

mod3.1 <- lm(likelinc~ one_noc + one_llc + one_rlc_1 + one_rlc_2 + one_rlc_3 + female + as.numeric(Q2.1) + edu + as.numeric(Q13_1) + income, df)
summary(mod3.1)

stargazer(mod3.1)

# Extract the coefficient and standard error for no coalition

coef3.1<-summary(mod3.1)$coefficients[,1] # coef
se3.1<-summary(mod3.1)$coefficients[,2] # se

mod.df3.1<-data.frame(cbind(coef3.1, se3.1))
mod.df3.1<-mod.df3.1[-c(1,3,4,5,6,7,8,9,10,11),]
varnames3.1<-c("No Coalition")

mod.df3.1$upper90 <- mod.df3.1$coef3.1 + qnorm(0.95)*mod.df3.1$se3.1
mod.df3.1$lower90 <- mod.df3.1$coef3.1 - qnorm(0.95)*mod.df3.1$se3.1
mod.df3.1$upper95 <- mod.df3.1$coef3.1 + qnorm(0.975)*mod.df3.1$se3.1
mod.df3.1$lower95 <- mod.df3.1$coef3.1 - qnorm(0.975)*mod.df3.1$se3.1
mod.df3.1$varnames3.1 <- varnames3.1

# Likelihood of formation of right-leaning coalition (LDP-Komeito-Nippon Ishin))

mod4.1 <- lm(likelilkn~ one_noc + one_llc + one_rlc_1 + one_rlc_2 + one_rlc_3 + female + as.numeric(Q2.1) + edu + as.numeric(Q13_1) + income, df)
summary(mod4.1)

stargazer(mod4.1)

# Extract the coefficient and standard error for right-leaning coalition

coef4.1<-summary(mod4.1)$coefficients[,1] # coef
se4.1<-summary(mod4.1)$coefficients[,2] # se

mod.df4.1<-data.frame(cbind(coef4.1, se4.1))
mod.df4.1<-mod.df4.1[-c(1,2,3,5,6,7,8,9,10,11),]
varnames4.1<-c("Right")

mod.df4.1$upper90 <- mod.df4.1$coef4.1 + qnorm(0.95)*mod.df4.1$se4.1
mod.df4.1$lower90 <- mod.df4.1$coef4.1 - qnorm(0.95)*mod.df4.1$se4.1
mod.df4.1$upper95 <- mod.df4.1$coef4.1 + qnorm(0.975)*mod.df4.1$se4.1
mod.df4.1$lower95 <- mod.df4.1$coef4.1 - qnorm(0.975)*mod.df4.1$se4.1
mod.df4.1$varnames4.1 <- varnames4.1

# Likelihood of left-leaning coalition (Komeito-DPFP-CDPJ-SDPJ-JCP

mod6.1 <- lm(likelikdcsj~ one_noc + one_llc + one_rlc_1 + one_rlc_2 + one_rlc_3 + female + as.numeric(Q2.1) + edu + as.numeric(Q13_1) + income, df)
summary(mod6.1)

stargazer(mod6.1)

# Extract the coefficient and standard error for left-leaning coalition

coef6.1<-summary(mod6.1)$coefficients[,1] # coef
se6.1<-summary(mod6.1)$coefficients[,2] # se

mod.df6.1<-data.frame(cbind(coef6.1, se6.1))
mod.df6.1<-mod.df6.1[-c(1,2,4,5,6,7,8,9,10,11),]
varnames6.1<-c("Left")

mod.df6.1$upper90 <- mod.df6.1$coef6.1 + qnorm(0.95)*mod.df6.1$se6.1
mod.df6.1$lower90 <- mod.df6.1$coef6.1 - qnorm(0.95)*mod.df6.1$se6.1
mod.df6.1$upper95 <- mod.df6.1$coef6.1 + qnorm(0.975)*mod.df6.1$se6.1
mod.df6.1$lower95 <- mod.df6.1$coef6.1 - qnorm(0.975)*mod.df6.1$se6.1
mod.df6.1$varnames6.1 <- varnames6.1

# Change names of different variables

names(mod.df4.1)<-names(mod.df3.1)
names(mod.df6.1)<-names(mod.df3.1)

# Merge the coefficients and standard errors from different coalition scenarios

mod.dfX1<-rbind(mod.df3.1, mod.df4.1, mod.df6.1)

varnames<-c("No Coalition","Right","Left")

# Plot

pX1<-ggplot() + 
  geom_linerange(data=mod.dfX1, mapping=aes(x=varnames, ymin=lower90, ymax=upper90), size=1.5) + 
  geom_linerange(data=mod.dfX1, mapping=aes(x=varnames, ymin=lower95, ymax=upper95), size=0.8) + 
  geom_point(data=mod.dfX1, mapping=aes(x=varnames, y=coef3.1), size=3) + 
  geom_hline(yintercept=0, colour="red", size=0.5) + 
  labs(x="Variable Names", y="Coefficient", title="Coalition Likelihood: Komeito") +  # Labels
  coord_flip() +  # Rotate the plot
  theme_bw()  # Nicer theme
pX1

### Experiment 2: DPFP

## Figure 4

# Likelihood of no coalition formation

mod8.1 <- lm(likelinc~ two_noc + two_llc + two_rlc + female + as.numeric(Q2.2) + edu + as.numeric(Q13_1) + income, df2)
summary(mod8.1)

stargazer(mod8.1)

# Extract the coefficient and standard error for no coalition

coef8.1<-summary(mod8.1)$coefficients[,1] # coef
se8.1<-summary(mod8.1)$coefficients[,2] # se

mod.df8.1<-data.frame(cbind(coef8.1, se8.1))
mod.df8.1<-mod.df8.1[-c(1,3,4,5,6,7,8,9),]
varnames8.1<-c("No Coalition")

mod.df8.1$upper90 <- mod.df8.1$coef8.1 + qnorm(0.95)*mod.df8.1$se8.1
mod.df8.1$lower90 <- mod.df8.1$coef8.1 - qnorm(0.95)*mod.df8.1$se8.1
mod.df8.1$upper95 <- mod.df8.1$coef8.1 + qnorm(0.975)*mod.df8.1$se8.1
mod.df8.1$lower95 <- mod.df8.1$coef8.1 - qnorm(0.975)*mod.df8.1$se8.1
mod.df8.1$varnames8.1 <- varnames8.1

# Likelihood of formation of right-leaning coalition (DPFP-Komeito-Nippon Ishin))

mod10.1 <- lm(likelidkn~ two_noc + two_llc + two_rlc + female + as.numeric(Q2.2) + edu + as.numeric(Q13_1) + income, df2)
summary(mod10.1)

stargazer(mod10.1)

# Extract the coefficient and standard error for right-leaning coalition

coef10.1<-summary(mod10.1)$coefficients[,1] # coef
se10.1<-summary(mod10.1)$coefficients[,2] # se

mod.df10.1<-data.frame(cbind(coef10.1, se10.1))
mod.df10.1<-mod.df10.1[-c(1,2,3,5,6,7,8,9),]
varnames10.1<-c("Right")

mod.df10.1$upper90 <- mod.df10.1$coef10.1 + qnorm(0.95)*mod.df10.1$se10.1
mod.df10.1$lower90 <- mod.df10.1$coef10.1 - qnorm(0.95)*mod.df10.1$se10.1
mod.df10.1$upper95 <- mod.df10.1$coef10.1 + qnorm(0.975)*mod.df10.1$se10.1
mod.df10.1$lower95 <- mod.df10.1$coef10.1 - qnorm(0.975)*mod.df10.1$se10.1
mod.df10.1$varnames10.1 <- varnames10.1

# Likelihood of left-leaning coalition (DPFP-CDPJ-SDPJ-JCP)

mod12.1 <- lm(likelidcsj~ two_noc + two_llc + two_rlc + female + as.numeric(Q2.2) + edu + as.numeric(Q13_1) + income, df2)
summary(mod12.1)

stargazer(mod12.1)

# Extract the coefficient and standard error for left-leaning coalition

coef12.1<-summary(mod12.1)$coefficients[,1] # coef
se12.1<-summary(mod12.1)$coefficients[,2] # se

mod.df12.1<-data.frame(cbind(coef12.1, se12.1))
mod.df12.1<-mod.df12.1[-c(1,2,4,5,6,7,8,9),]
varnames12.1<-c("Left")

mod.df12.1$upper90 <- mod.df12.1$coef12.1 + qnorm(0.95)*mod.df12.1$se12.1
mod.df12.1$lower90 <- mod.df12.1$coef12.1 - qnorm(0.95)*mod.df12.1$se12.1
mod.df12.1$upper95 <- mod.df12.1$coef12.1 + qnorm(0.975)*mod.df12.1$se12.1
mod.df12.1$lower95 <- mod.df12.1$coef12.1 - qnorm(0.975)*mod.df12.1$se12.1
mod.df12.1$varnames12.1 <- varnames12.1

p12.1<-ggplot() + 
  geom_linerange(data=mod.df12.1, mapping=aes(x=varnames12.1, ymin=lower90, ymax=upper90), size=1.5) + 
  geom_linerange(data=mod.df12.1, mapping=aes(x=varnames12.1, ymin=lower95, ymax=upper95), size=0.8) + 
  geom_point(data=mod.df12.1, mapping=aes(x=varnames12.1, y=coef12.1), size=3) + 
  geom_hline(yintercept=0, colour="red", size=0.5) + 
  labs(x="Variable Names", y="Coefficient", title="Coalition Likelihood: DCSJ Coalition (w/ covariates)") +  # Labels
  coord_flip() +  # Rotate the plot
  theme_bw()  # Nicer theme
p12.1

# Change names of different variables

names(mod.df10.1)<-names(mod.df8.1)
names(mod.df12.1)<-names(mod.df8.1)

# Merge the coefficients and standard errors from different coalition scenarios

mod.dfX2<-rbind(mod.df8.1, mod.df10.1, mod.df12.1)

varnames<-c("No Coalition","Right","Left")

# Plot

pX2<-ggplot() + 
  geom_linerange(data=mod.dfX2, mapping=aes(x=varnames, ymin=lower90, ymax=upper90), size=1.5) + 
  geom_linerange(data=mod.dfX2, mapping=aes(x=varnames, ymin=lower95, ymax=upper95), size=0.8) + 
  geom_point(data=mod.dfX2, mapping=aes(x=varnames, y=coef8.1), size=3) + 
  geom_hline(yintercept=0, colour="red", size=0.5) + 
  labs(x="Variable Names", y="Coefficient", title="Coalition Likelihood: DPFP") +  # Labels
  coord_flip() +  # Rotate the plot
  theme_bw()  # Nicer theme
pX2
